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ABSTRACT 


The goal of this study was to experimentally evaluate 
and mathematically model the performance of phase change 
thermal control devices containing high thermal-conductivity 
metal matrices. Three aluminum honeycomb fillers were 
evaluated at five different heat flux levels using n-oct- 
adecane as the test material. Initially phase change per- 
formance was evaluated with no filler in the n-octadecane 
so that a base line performance for each heat flux level 
could be established. 

The experimental equipment consisted of a test cell, 
two electric heaters, a watt meter, two ammeters, and a 
multipoint recorder. The test chamber measured 15.24-by- 
7.62-by-2.54 centimeters (6-by-3-by- 1-inches) . The cell 
was heated by two 7 * 62 -by- 7.62 centimeters ( 3-by-3-inch) 
electric heaters, which were held at a constant heat flux. 
The amp and watt meters provided the measurement of the 
heat flux to the heaters. Temperature responses to the 
upset were measured by 16 copper-constantan thermocouples 
located throughout the test cell. 

The system was mathematically modeled by approximating 
the partial differential equations with a three-dimensional- 
implicit-alternating direction technique. This implicit 
method was used so that the small time step required by 
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the presence of the metallic filler in the explicit method 
could be eliminated. The boundary conditions used in this 
model were a temperature profile on the copper heating plate 
and insulated on the other five sides. 

The mathematical model predicts the system quite well. 
All of the phase change times are predicted. The heating 
of the solid phase is predicted exactly while there is some 
variation between theoretical and experimental results in 
the liquid phase. This variation in the liquid phase could 
be accounted for by the fact that there are some heat losses 
in the cell and there could be some convection in the ex- 
perimental system. 
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INTRODUCTION 


Phase-change phenomena have received wide scientific 
attention for some time. Phase change is of significant 
importance in many technical problems such as solidifi- 
cation of an asphalt layer, melting and solidification of 
metals and alloys and general crystal growth. In recent 
years solid-liquid phase change has been used for thermal 
control devices in space vehicles. In concept, such 
materials would be used in passive systems that employ 
the process of melting or solidification to remove or add 
thermal energy from or to a system. 

In a study by Northrop Corporation ( 1 ) , the prop- 
erties that phase-change materials must have in order to 
control the temperature of electronic equipment were found. 
The phase-change material should be nontoxic, chemically 
inert, stable and noncorrosive. The material should also 
have small density variations with a high latent heat of 
fusion. The material should also melt in 283 to 338°K 
(50-to 150°F) range. N-paraffins with an- even number of 
carbon atoms are the most promising phase-change materials. 
N-octadecane was used in this study. 

Virtually all of the currently used phase-change 
materials have a low thermal diffusivity. Therefore, 
their use in phase-change thermal control units is hampered 
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by their inability to conduct heat. Their thermal con- 
ductivities are on the same order of magnitude as thermal 
conductivities of some of the best insulating materials. 

A method to improve heat transfer rate of the phase-change 
material is to surround the phase-change material with a 
high thermally conductive metal matrix. The metal would 
conduct the heat and the phase-change material would 
absorb the heat load. 

Bentilla, Sterrett and Karre (l) evaluated a number 
of metallic fillers in a phase-change environment. They 
evaluated aluminum foam, aluminum wool, aluminum honeycomb 
and copper foam. Their work showed that the most advanta- 
geous type of filler was the aluminum honeycomb. Hale, 
Hoover and O'Neill (2) developed a method using an overall 
energy balance to predict the performance of thermal control 
phase-change devices. In their study they found that by 
neglecting three-dimensional heat transfer significant 
errors were introduced. The goal of this study was to 
increase the understanding of increased thermal diffusivity 
phase-change devices by experimentally evaluating and 
mathematically modeling them. With the aid of the three- 
dimensional mathematical model developed in this study 
engineers will be able to design phase-change thermal 
control units more efficiently. 



LITERATURE SURVEY 


There has been a large amount of literature pub- 
lished on the subject of phase-change phenomena . This 
literature survey deals with only a small portion of the 
published material. One of the first studies of phase 
change was made by Carslaw and Jaeger (3)* In their study 
they developed an approximate mathematical model for semi- 
infinite bodies. They discussed the problem of modeling 
a system with a moving interface. No exact solutions 
were given for the mathematical modeling of finite bodies 
with phase change. 

At this institution two studies have been completed 
which are concerned with the one-dimensional interface 
equation given by Arpaci (4). In the first study Pujado, 
Stermole and Golden (5) developed a theoretical model for 
the one-dimensional melting of a finite paraffin slab. 
Finite difference methods were used to solve the partial 
differential equations governing the physical system. 

This model solved the two-phase one-dimensional heat- 
transfer equations with phase change and variable thermal 
properties. In their theoretical analysis Pujado, Stermole 
and Golden neglected free convection in the liquid phase. 
The results from the study were in agreement with those of 
an earlier study by Northrop Corporation. The second study 
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by Ukanawa, Stermole and Golden (6) was the investigation 
of the solidification of a finite amount of liquid para- 
ffin. A mathematical model was developed to solve the 
two-phase heat-transfer equation with phase change. This 
model used constant thermal properties for each phase and 
moveable boundary conditions. The model neglected convec- 
tion, supercooling and nucleation effects. The comparison 
was good between theoretical and experimental results. 

There have been three studies at this institution 
which concern the two-dimensional phase-change problem. 

Shah (7) investigated the solidification of r.-octadecane 
using microphotographic equipment and a temperature re- 
corder. In this study a two-dimensional mathematical model 
was developed to predict the temperature profiles in the 
freezing paraffin and the average interfacial height during 
the solidification process. The model neglected convection 
An approximate method was used to calculate the phase 
change. A presentation of the various types of phase- 
change calculations is given by Dusinberre (8). Reasonably 
good agreement was obtained between the experimental and 
theoretical results. The second study by Bain, Stermole 
and Golden (Q) was an investigation of the gravity-induced 
free convection in the melting of a finite paraffin slab. 
Temperature profiles were measured when the test cell was 
inclined at different angles to produce the free convection 
In this study a two-dimensional pure-conduction model was 
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presented. Finite difference methods were used to solve 
the governing partial differential equations while the 
method of excess degrees was used to calculate the phase 
change. Since the model neglected convection there were 
significant differences between the theoretical and 
experimental results. This study showed that gravity- 
induced free convection can be important in the melting 
process. The third investigation by Ukanawa ( 10) studied 
the effect of gravity- induced free convection upon the 
solidification of a finite paraffin slab. A two-dimen- 
sional heat-transfer model was developed in this study. 

The model used an important velocity profile and a 
limiting velocity in the convection calculation. A pseudo 
heat-capacity was used to calculate the phase change. 

Other papers have also been published which deal with 
the melting of finite slabs. Chi-Tien and Yin-Chao Ten (ll) 
presented an approximate solution for the temperature dis- 
tributions and melting rate. The heat transfer was by 
natural convection caused by buoyancy forces. Goodman and 
Shea (12) developed a series solution to solve the one- 
dimensional problem of the melting of a finite slab. 

Crank and Nicholson (13), Douglas and Rachford ( 1 4 ) , 
Peaceman and Rachford (15) -md Brian ( 16) all developed 
three-dimensional finite difference techniques to solve 
the unsteady state heat-conduction partial differential 
equation. The three-dimensional alternating direction 
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technique developed by Brian was used in this study. 

Brian showed this technique is unconditionally stable and 
has the accuracy of the Cranlc-Nicholson method. 

Grodzka (17) suggested two methods to increase the 
thermal dif f usivities of phase-change materials and thus 
the heat-transfer rate into the thermal control unit. One 
ivas to put metallic panels in with the phase-change material. 
The other was to mix the phase-change material with a com- 
patable but higher thermal diffusivity material. She pursued 
the second method and suggested several possible materials 
that could be mixed with the phase-change material. In 
their recommendation for the improvements of thermal control 
phase-change packages, Shlosinger and Bentilla (18) 
recommended that a metallic filler be added to increase 
the heat transfer rate. Bentilla, Sterrett and Karre (!) 
evaluated a number of metallic fillers in a phase-change 
environment. These fillers included aluminum wool, 
aluminum foam, aluminum honeycomb and copper foam. Their 
work showed that the aluminum honeycomb was the most 
advantageous geometry for a high thermal diffusivity 
filler material. Hale, Hoover and O'Neill (2) developed 
a method by which the ratio of the phase-change material 
to the filler material may be optimized. They used an over- 
all energy balance to optimize the ratio of the filler 
material area to phase-change material area as a function 
of the hot plate temperature. 



THEORY 


The finite difference equation which can be used to 
mathematically model a nonhoinogeneous system with phase 
change will be developed in this section. The test cell 
was heated from the top to minimize convection, however, 
presence of the metal matrix is likely to cause convection. 
The paraffin closer to the metal matrix will heat up 
faster due to the high rate of heat transfer through the 
metal and this situation will cause some convection. The 
convection caused in this way will be considered negligible. 
Since the assumption was made that convection currents 
have no significant effect on the heat transfer in the 
test chamber, this development will neglect convection. 

For a discussion of mathematically modeling a phase-change 
system with convection see reference 10. The general three- 
dimensional heat-conduction model will be developed first, 
then the three-dimensional alternating direction technique 
of Brian (l6) will be discussed. Finally, the general 
heat-conduction equation and the three-dimensional alter- 
nating direction technique will be applied to the non- 
hoinogeneous phase-change problem in this study . 

General Three-Dimensional Equation 

The general heat-conduction equation may be derived 
by making an energy balance on a three-dimensional non- 
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homogeneous nodal system. Making an energy balance on 
node i,j,k in figure 1 yields 

Energy in - energy out + energy generated 
= energy accumulated (l) 


or , 


( P C p V) eff (T*(i,j,k)-T(i,j,k)) = (q x A x ) irj -(q x A x ) out 


J 1 

At 


+ ( q A ) . - ( q A ) , 

y y xn y y out 


+(qA ). - ( q A ) , 

M z z in z z out 


+ GV 1 


( 2 ) 


where 




A 


z 


G 

V' 


T* 


T 

t 


= the heat flux in the x-direction 
= the heat flux in the y-direction 
= the heat flux in the z-direction 
= the cross sectional area of the node i,j,k 
perpendicular to the x-direction 
= the cross sectional area of the node i,j,k 
perpendicular to the y-direction 
= the cross sectional area of the node i,j,k 
perpendicular to the z-direction 
= energy generated per volume 
= volume of the material generating energy 
= the temperature at time t+At 
= the temperature at time t 
= time 




Figure 1 


Three-Dimensional Koclal System 
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At = time increment 

(^ C pV) e ff = Che e f^ ec ^^- ve (pC pV) for a nonhomo- 
geneous system 

The heat flux may be found in terms of the temperature 
difference of the nodal system from Fourier's law. From 
Fourier's law one obtains 

r ■ 

(q x ) in = K(T(i-l,j,k)-T(i,j,k))/Ax (3) 

where 

K = the thermal conductivity of the material 
Ax = the incremental distance between nodes in the 
x-direction 

Writing a similar equation for each of the other heat flux 
terms in equation (l) and substituting back into equation 
(l) yields 

( P c p v) eff (T*(i,j,k)-T(i,j,k)) = 

A t 

(( KA ) x ) e ff (T(i-1, j,k)-T(i,j,k) )/Ax - 
((KA) x ) eff (T(i, j,k)-T(i+l, j,k) )/Ax + 

((KA) y ) eff (T(i, j-l,k)-T(i,j,k) )/Ay - 
( ( KA ) y ) e f f (T(i, j,k)-T(i, j+l,k ) ) / AY + 

((KA) z )eff (T(i, j,k-l)-T(i, j,k) )/Az - 

( (KA) z ) eff (T(i,J^k)-T(i,j,k+l))/Az + GV' (4) 


' I 


. .'I i i 


((KA) x ) e ££ = the effective thermal conductivity times 
the effective area in the x-direction for 
a nonhomogeneous system 

((KA) ) = the effective thermal conductivity times 

v y ef f 

the effective area in the y-direction for 
a nonhomogeneous system 

( ( KA ) ) ~~ = the effective thermal conductivity times 

" .z err 

the effective area in the z-direCtion for 
a nonhomogeneous system 

Ay = the incremental distance in the y-direction 

A 7 == the incremental distance in the z-direction 

Equation (4) is the general three-dimensional heat- 

conduction equation written in an explicit finite difference 

form. To use this equation one must define and evaluate the 

effective PC V and the effective KA terms. The effective 
r P 

pc V will be defined as the sum of the p C V terms repre- 


senting each material present in the node or, 

N 

(pc v) = Z P C V 

'* p 'eff 'n pn n 

where 


( 5 ) 


N = the number of different materials present in the 
node 

This is just the sum of the heat capacitance of each of 
the different materials. An expression for the effective 
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KA term can be found by looking at the one-dimensional 
steady state energy balance on a nonhomogeneous node. A 
three-dimensional system could be used to derive this 
expression, but the analysis is simpler in one-dimension 
while the same result is obtained. An energy balance on 
node m in figure 2 yields 
(qA) ln -(qA) out = o 


( 6 ) 


The heat flux into node m is the sum of the heat which flows 
through each of the different materials, or 

( qA ) AjQj + ^ 2^2 ^ 3^3 ( 7 ) 


where 

— the heat flux through material 1 
= the heat flux through material 2 

q^ = the heat flux through material 3 

Aj = the area of material 1 perpendicular to the 
heat flux 

A 2 = the area of material 2 perpendicular to the 
heat flux 

A^ = the area of material 3 perpendicular to the 
heat flux 

From Fourier's law the heat flux terms may be evaluated as 
f ollows : 

( q A > in = < KA) eff (T l~ T m )/Ax (8) 

* 

A 1 q 1 = ( KA )i ^ T 1 1 _T lm^A x 


( 9 ) 
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Figure 2 

One-Dimensional Fonhomogeneous • 
Kodal System 
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A 2 q 2 ( KA )2 21~ T 2w) ^ AX 

A 3 q 3 = (KA) 3 (T 31- T 3m } / AX 


( 10 ) 

(ID 


Substituting equations (8) through (ll) into equation (7) 
yields 


(KA) eff (T r T m ) = (KA) i (T 11 "T lni ) + (KA) 2 (T 21 -T 2m ) 

+ (KA) 3 (T 31 -T 3m ) (12) 

Since there is assumed to be no gradients in a node in a 
finite difference network the following is true. 



T i 

= T 11 

T 21 

" T 31 


T 

m 

= T lm 

T 2m 

T 3«i 


The effective KA may be found by substituting equations 
and (14) into equation (12) and dividing by T^-T^ This 
procedure yields 


(13) 

(14) 
(13) 


(KA) eff = (KA) 1 + (KA) 2 + (KA) 3 (15) 

or more generally 

(KA) eff = ZK n A n (16) 

For the nonhomogeneous system we must define an effective 
pC V for each node and an effective KA for each of the six 
sides for each node. Equation (4) may now be rewritten in 
terms of the above relations as 


T*(i, j ,k)-T(i, j ,1c) = 

At(KA(i,.i,k) x ) eff (T(i-1, j,k)-T(i, j,k) ) 
4x(/3 c p v (i, j/k) ) eff 
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- At(KA(i+l, j,K') x ) eff (T(i, j,k)-T(i+l, j,k) ) 

Ax(pCpV(i, j,k) ) eff 

+ At(KA(i, j,k) y ) eff (T(i,j-l,k)-T(i, j,k) ) 

Ay(p c p v (i, j,k) ) eff 

- At(KA(i, j+l,k) y ) eff (T(i, j,k)-T(i, j+l,k)) 

Ay(^)C p V(i, j,k) ) eff 

+ At(KA(i, j,k) z ) eff (T(i, j,k-l)-T(i,j,k) ) 

^z(p C p V ( i ,J, k )) e ff 

- At(KA(i, 3 ,k+l) z ) eff (T(i, j,k)-T(i, j,lc+l)) 

Az(pCpV(i, j,k) ) eff 

+ GV 1 At 

C p V(i, j } k) eff (17) 


where 


(KA ( i f 3 } k) x ) 


(KA( i+1 j3j k ) x ) e ££ 


(KA( i , 3 j k ) y ) e f f 


(KA(i,j+l,k) y ) eff = 


the effective KA term perpen- 
dicular to the x-direction on the 
left side of the node 
the effective KA term perpen- 
dicular to the x-direction on the 
right side of the node 
the effective KA term perpen- 
dicular to the y-direction on the 
left side of the node 
the effective KA term perpen- 
dicular to the y-direction on the 
right side of the node 
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(KA( i, j , k) z ) ^ = the effective KA term perpen- 

dicular to the z-direction on the 
left side of the node 

(KA( i, j , k+1 ) z ) = the effective KA term perpen- 
dicular to the z-direction on the 
right side of the node 

pCpV( i, j , k) = the effective pC p V term for the 

node 

A simpler form of equation (17) can be written in terms of 
variable constants in place of the effective terms. 

T- ;: (i,d,k)-T(i } 3,k) = Cj j (T( i- 1 , j , k) -T(i, j , k) ) A t 

- C 12 (T(i,a,k)-T(i+l,j,k))At 

+ C 21 (T(i,j+l,k)-T(i,j,k))At 

- C 22 (T(i,j,k)-T(i,j+l,k))At 

+ j ( T ( i, 3 , lc- 1 )-T( i, j , k) ) At 

- C 32 (T(i,j,k)-T(i,j,k+l))At + GC 4 At ( 18 ) 

where the constants are defined as follows ar.d are evaluated 


for each node as indicated. 

C j j ~ (KA( i ? j f k)^) e -p-f 

- 2 KA(i,j,k,n) x 


P ^"p^ ( -*-» J f k) Ax 

4xIpC p V(i,j,k,n) 

(19) 

C J2 — ( KA ( i+ 1 j j , k ) x ) 

= ZKA(i+l,.j,k,n) x 


P c p V(i,.l,k) eft dx 

A''Z/>C p V( i, 3 ,k,n) 

(20) 

C 21 = (KA(i.j,k) y ) eff 

= ZKA(i,j ,k,n) y 


P f pV ( i > 3 y k ) e p f Ay 

W LP C p v ( i »J»k,n) 

(21) 
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C 22 = (KA(i,j+l,k) y ) eff 

= Z KA(i, j+r,k,n) y 


/>CpV(i, j,k) eff Ay 

zp c p v ( ± > j > k > n ) 

(22) 

C 31 = ( KA ( x , j , k ) ^ ) e f f 

= Z KA(i, j ,k,n ) ? 


pC p V(i, j,k) eff A.z 

Z PCpV (i,j,k,n) 

(23) 

C, 2 = (KA(i,j,k+I) z ) et . f 

= Z K A ( i , j , k+ 1 , n ) z 


pC p V(i, j,k) eff Az 

Az JC/lCpV (i,j,k,n) 

(24) 

c 4 = V'/(pC p V(i,j,k) eff 

= V * 



p c p V(i, j,k,n) 

(25) 

Equation (l8) may be used 

over the entire nodal 

system 


by just varying the above constants in order to define the 
effective terms for each node. By redefining the constants 
in equation (l8) it may be used on the boundaries also. To 
see how this may be done look at the one-dimensional form 
of equation ( 1 8 ) . Again a three-dimensional analysis will 
give the same result, but the one-dimensional analysis is 
simpler . 

T*(i, j,k)-T(i, j,k) = C n (T(i-1, j,k)-T(i, j,k) ) At 

+ C l2 (T(i, j,k)-T(i+l, j,k) ) At 

+ GC^At (26) 

Figure 3 shows the four types of boundary conditions 
that can be used. Each type will be described in the 
following section. 

Type I : This boundary condition represents a system 

in which there is heat transfer to a fluid. An energy 
balance on node 2 of this nodal system yields 
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T* 2 -T 2 = At hA(Tj-T 2 ) 

- At KA(2) eff (T 9 -T 3 ) 


P° p V < 1) eff 

**P C P v( ')eff 

(27) 


Comparison of equation (27) to equation (26) shows that 
C 1 = hA 

P C p V<1, eff ( 28 ) 

C . = 0; since there is no generation of ( 29 ) 

energy in this node 

Cj 0 is defined by equation (20). 

Typos TT and III : These are two types of insulated 

boundary conditions. An energy balance on node 2 of these 


two systems yields 

T* -T, ; = -At KA(2) ff (T 9 -T ) 

^V (1, eff ’ (30) 

Comparing equation (30) with equation (26) indicates that 

C U ° (31) 

G = 0 (32) 


T} r pe TV : This type of boundary is insulated and 

energy is being generated in the boundary node. An energy 
balance on this system yields 


T *2~ T 2 = _At KA(2) eff (T 2“ T 3 ) + tJ " V ' Afc 

4x / JC p V(1J eff P C p V < 1, eff (33) 

A comparison of this equation with equation (26) yields 

C n = 0 (34) 

G = U" (35) 
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By setting G equal to the energy source in any node in 
the system, there may be a generation term in any node 
or boundary equation. 

Th ree- Dime n sional Alternating Direction Technique 

The time step required by the stability criteria in 
the explicit solution of equation ( 1 8 ) is extremely small 
due to the presence of the metallic matrix. The required 
time step is such that the computer time required to model 
an experimental run exceeds the actual experimental run 
time. Therefore, an implicit technique was used to elim- 
inate the stability requirement on the size of the time 
step. The implicit method used in this study is the 
three-dimensional alternating direction technique developed 
by Brian (l6). This method is a variation of the Douglas- 
Rachford method that has the accuracy of the Crank- 
Nicholson method and has been shown to be unconditionally 
stable . 

The three-dimensional technique of Brian solves for 
three intermediate half-time-step temperatures and then 
uses these to solve for the new f ull-time-step temperature. 
The equations that demonstrate this method are as follows: 

A^(CT*(i,j,k)) +A^(CT(i,3,k)) + A^(CT(i,j,k)) = 

X j ^ 


T*-(i,.i . k ) - T ( i , ,j ,k) 

At/ 2 


( 36 ) 
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where 


A 2 (CT*(i,j,k)) + A 2 (CT**(i,j,k)) + A 2 (CT(i,3,k)) = 
x y 

T**(i..i ,k) 

At/2 (37) 


A~. ( CT* ( i , j , k ) ) + A^CT**(i,j,k)) + A ( CT*** ( i , j , lc ) ) = 

X V . 

(38) 


T**»(l, j,k)-T(i,j,k) 
At/2 


A * ( CT* ( i , j , k ) ) + A* ( CT** ( i , j , k ) ) +A l z ( CT*** ( i , 3 , k ) ) 


T 1 (i,;i , k ) - T ( i , ,j ,k) 
A t 


(39) 


A 2 (CT(i,j,k)) = C n T(i-l, j,k)~(C 11 +C 12 )T(i, j,k) 

+ C 12 T(i+l,j,k) (40) 

A 2 (CT(i, j,k)) = C 21 T(i, j-l,k)-(C 21 +C 22 )T(i, j,k) 

+ C 22 T(i,j+l,k) (41) 

A 2 (CT(i,j,k)) = C 31 T(i,.i,k-l)-(C 31 +C 32 )T(.i,.i,k) 

+ C 32 T(i, j,k+l) ( 42 ) 

T = the temperature at time t 

T' = the temperature at the new full time step t+At 
T* , T** , T*** = the intermediate half-time-step 
temperatures 
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C,,, C to , C ', C 01 and C 00 are evaluated for each node 
and are defined by equations (19) through (24). For 
stability it is imperative that T* be used in the x- 
direction difference, T** be used in the y-direction 
difference and T**# be used in the z-direction difference. 

A simpler form of the equations (36) through (39) 
were used to solve the problem presented in this thesis. 

When T*-*-” is eliminated from equations (36) through (39) , 
the following set of equations results . 

A^(CT*(i>j,k)) + Ay(CT(i,j,k)) + A^(CT(i,j,k)) = 

T*(i, j 

At/2 (43) 


Ay(CT**(i,j,k))- Ay(CT(i,j,l<)) = 

T* * ( i , ,j ,k) -T#(i, j , k) 

At? 2 


(44) 


^(CT«(i,j,k))- A^(CT(i,j,k)) = 

(45) 

Equation (44) is the difference between equations (36) 
and (37), while equation (45) is found by eliminating 
T-”-** from the difference of equations (37) and (38) and 
the difference of equations (38) and (39). Equation (43) 
relates the unknown T# values along a row parallel to the 
x-axis. Y/hen equation (43) is solved for the unknown T* 



values a system of simultaneous equations results: 
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b l T l* + 

O 

h-k 

H 

to 


il 

d l 

a 2 T l* + 

b 2 T 2 r + C 2 T 3 ' 


= 

d 2 


a 3 T 2 * + b 3 T 3 * + 

c T * 

3 4 

= 

d 3 

a 0 T* 

n-2 n- 

-3 + b n-2 T ”n- 2 + 

c 0 T* . 

n-2 n- 1 


d n- 2 


a n-l T *n-2 + 

b , T'“ , 

n- 1 n- 1 

= 

d n- 1 


( 46 ) 


The values of the coefficient and c^ are the 

coefficients of the unknown T* temperatures and d. is 

the sum of the remaining terms. It is supposed that the 

I 

grid points are' designated by 0, 1 , 2 , 3 , . • • n- 1 , n and 

that and T are determined from the boundary conditions. 
0 n 

The matrix of the coefficients a, b, and c is 
tridiagonal. There is a very efficient method of the 
solution for the tridiagonal system. The value of T^* 
in equation (46) can be found by following procedure: 


T* . = F . 
n- 1 n- 1 


T*. = F.-c.T* M /w. 

x xx n+r x 


( 47 ) 

( 48 ) 


where W. and F. are determined by the following recursion 
xx 


formula 


w. = b.-a.c. 1 with w. = b. 

x x 1 x- 1 1 1 


IV 


i- 1 


F . = d . -a . F . 1 

x x x x- 1 


with F^ = dj/bj 


( 49 ) 


( 50 ) 
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Equation (44) is then solved in a similar manner, 
but this time the simultaneous equations relate T** values 
along a row parallel to the y-axis . The solution of 
equation (45) is then found in a similar manner with sets 
of tridiagonal simultaneous equations relating T***'* 
values along a row parallel to the z-axis . 

Nonhomogeneous Phase-Change Problem 

The system in this study consists of a hexagonal 
aluminum matrix in n-octadecane . The physical prop- 
erties of aluminum and n-octadecane are given in ref- 
erence 19 and 20 and are tabulated in table 1. Due to 
the symmetry of the system only a small portion of the 
test cell must be modelled. Figure 4a shows how the 
test cell can be broken down by lines of symmetry in the 
two horizontal directions. Using these lines of symmetry 
the filler system can be broken down into the system 
shown in figure 4b. The sides in the x and y directions 
are considered insulated from lines of symmetry. The 
bottom of the cell will be considered insulated in this 
analysis. Actually this may not be the case but fairly 
good agreement between theoretical and experimental data 
results if the bottom is considered insulated in the 
theoretical analysis. The nodal system will be defined 
as shown in figures 5a and 5b. The x-direction is indi- 
cated by i, j indicates the y-direction and k represents 
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Table 1 

Literature values of the physical properties for n-oct- 
adecane and aluminum. 

N-OCTADECAKS 

Density 

p (solid) = (-0.0003336)T + 1.0918; g/cc 

P ( liquid) = (- 0 . 00 1 2505 ) T + I.I 3165 g/cc 

Heat Capacity 

C p (solid) = 2.164; watt-sec/ gni/°K 

C p (liquid) = (O.OOS213)T - 0.142375 watt-sec/gm/°K 

Conductivity 

IC (solid) = (-0.50054 x 10“ 5 )T + 0.002914; watt/ cm/°IC 

IC (liquid) - (-0.50054 x 10 _5 )T + 0.002914; watt/ cm/°K 

Melt point = 300.60 °K 

H^ (liquefaction enthalpy) = 243»9; watt-sec/ gm 
ALUMINUM 

p = 2.685 giu/cc 

Cp = 0.9792 watt-sec/gm/°K 

K = 0.1282 watt/cm/°K 



Figure 4a 


Aluminium Honey Comb 
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the z— direction . When the system of equations (43) 
through (44) j s applied to equation (l8) the following 
systems of equations results, where At = At/2 

For all i., j and k the a.^, b_^, c_^ and d. are 
defined as follows for equation (43)* 


a 


i 



At 


(5D 


b i = 1 + (C U + c i2 )4t 


(52) 


c i “ C 12 At 

d-i = T ( ± >^ k ) + (4y(CT(i,d,k)) +^(CT(i,j,k)))At 


(53) 

(54) 


^12* ^21 > C 3 1 anf l are calculated for each node 

according to equations ( 19 ) and (24). If there is no 


filler in the node the volume and. area of the filler in 
these calculations is set equal to zero. The insulated 
boundary conditions are evaluated as follows: 
when i and/or j = 2 for 3< k < K 

C U " 0 

(55) 


when 


C 21 = 0 

i = I and/or j = J for 3<;k <K 


(56) 



= 0 


when k — K for 2<;j<J and 2-C3‘.<;I 


C 


32 


0 


( 57 ) 
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A temperature boundary condition was used on the copper 
heating plate in this study. This condition is taken care 
of by setting T(i,j,2) equal to the heating plate temp- 
erature. 

For all i, j and k the a., b . , c. and d. are defined 

J J J J 

as follows for equation (44), 

a. = -C 91 At (58) 

3 21 


b. - 1 + (C 21 + C 22 )At 


c . = -C 0 

i 1Z 


(59) 

( 60 ) 
( 61 ) 


cl. *= T'”* ( i , j , k ) - (CT(i,j,k ))) At ( 6 l 

The boundary condition coefficients are calculated by 
equations (55) through (57). 

For all i, j and k the a k , b k , c k and d k are defined 

as follows for equation (45). 


a k ~ C 31 At 


b k = 1 + (C 31 + C 3 2 ) 4 1 


c k ~ C 32 A " 


( 62 ) 

(63) 

(64) 


d k = 2T**(i,j,k)-T(i, j,k)-(A^(CT(i, j,k)) + E)4t (65) 


where 


E =0 for 2<i<I; 2<,jCJ; 4<kC « 

E = C 3 j T ( i , ,j , 2 ) for 2<i<I; k=3 

c k = 0 for 2 <i<I; jc J ; k=3 


(67) 


( 68 ) 
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A summary of the various phase-change calculation 
techniques is given by Dusinberre (8). A modification 
of the method of excess degrees was used in this study. 

Since the n-octadecane used in this study was not a pure 
material (practical grade) it was allowed to melt over a 
1 . 76°K ( 3°F) temperature range. Since the heat capacity 
of n-octadccane is the same above and below the melting 
point, a term which has the units of degrees results 
when the latent heat is divided by the heat capacity. 

This term is called the excess degrees, which is the 
number of degrees the node would have risen if the phase 
change had not occurred. The following procedure is used 
to calculate the phase change. 

T(i,.j,k) .R.T m (69) 

If T(i,j,k) < T mo (70) 

the node is still solid. 

If T(i,j,k) > T mo (71) 

the node is in the process of melting. 

If equation (71) holds then the following procedure is 
followed. 

T e (i,.j,k) = (T(i,j,k)-T mo ).R.H f /C p (72) 

If T e (i,j,k) < H f /C p (73) 

the node is partially melted and its temperature 


is defined by 
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T(i,j,k) = T no + T e (i,j,k) *C p *1.76/ H f (74) 

But if 

T e (i,j,l<) > H f/ C p (75) 

the node has melted and the temperature of the now 
liquid node is given by 

T(i,j,k) = T mo (i,j,k) + 1.76 + ■ 

(T e (i,j,k)*C p - H f )/C p (76) 

Each new full-time-step temperature is corrected by the 
above phase-change calculation. 

A computer program was written using the above 
procedure to solve the nonhomogeneous heat-transfer 
problem with phase-change presented in the study. The 
results are shown and described in the discussion of 


results . 



EXPE RIMENTAL EQUIPMENT AND PROCEDURE 


A description of the experimental equipment and the 
experimental procedure are given in this section. 

Equipment 

The equipment used in this study can be seperated 
into three sections which are the test cell, the power input 
measuring system and the temperature recording system. 

These sections are discussed below. 

Test Cel l: The test cell (figures 6 and 7) con- 

sisted of a rectangular test chamber and a heating plate. 

The test chamber, 7 . 62-by- 1 5 . 24-by-2 . 54 centimeters (3-by- 
6-by- 1- inches ) was milled out of a block of plexiglass to 
minimize the sources of leaks. The heating plate was a 
10. l6-by-17.78-by-0.625 centimeter ( 4 -by- 7 ~by- 0 . 25 -inch) 
copper plate upon which two 7.62-by-7-62 centimeter (3-by- 
3-inch) electric heaters were epoxied. The electric heaters 
were obtained from Electrofilm Incorporated of North Holly- 
wood California. The cell was sealed by compressing the 
0-ring. During the run the cell was encased in approximately 
3.81 centimeters (1.5 inches) of styrofoam. 

Power Input Measuring System : This system, shown in 
3, consisted of a seven and one half ampere 
powerstat, a Hickok watt meter and two Welch A.C. am- 
meters. The powerstat provided a variable source of power 
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to the electric heaters. This power was measured by the 
watt meter. The output of the watt meter is a direct 
measure of the heat flux liberated by the electric heaters. 
The function of the ammeters was to insure that each 
heater received the same current and thus provide even 
heating . 

Temperature Recording System ; Sixteen copper- 
constatan thermocouples and a Bristol multipoint recorder 
comprised the temperature recording system. The recorder 
was able to record each point every two seconds with an 
accuracy of + 0.4267°K (+ 0.75°F). The thermocouple wires 
in the test cell were encased in glass probes. The glass 
probes served two purposes . One was to insulate the 
thermocouple wire from the metallic filler, the other was 
to keep the thermocouple at a constant height. The thermo- 
couple locations are given in table 2 , where coordinate 
0, 0, 0 is the left front upper corner of the test chamber . 

Exp e rlme n tal Procedure 

Experimental runs were made using the following 
procedure . 

1 . The cell Avas filled with n-octadecane after the 
filler material had been put into the test chamber. The 
cell was then bolted down to seal it. 

2. The cell was leveled to help minimize convection. 

3 . The recorder was turned on to record the initial 


P 
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temperature of the cell. 

4. When the test chamber had reached a constant 
initial temperature the powerstat was turned on to start 
the run . 

5 . The run was allowed to continue until the hot 
plate reached a temperature of 338.6l°K (150°F). At this 
temperature the plexiglass began to deform around the 
copper plate. 
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Figure 6 

Test Cell - Exploded View 
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Table 2 


THERMOCOUPLE LOCATIONS 

The location of the coordinate 0,0,0 is the left front 


upper corner 
the 15.24 cm 
is along the 
the vertical 

Thermocouple 

of the test 
side of the 
7. 62 cm side 
2.54 cm. 

Me . 

chamber. The x-axis is along 
test chamber, the y-direction 
while the z- direction is a Ion, 

x y z 

cm cm cm 

1 . . 


« • ♦ 

1.905 

.9525 

0.635 

2 . . 


• • • 

5.715 

.9525 . 

0.635 

3 • • 


t • * 

9.52 5 

-9525 

0.635 

4 • • 


• • • 

. 1 3-335 

.9525 

0.635 

5 • • 


• • • 

1.905 

2.8575 

1.27 

6 • • 


• • • 

5.7 15 

2.8575 

1.27 

7 • • 


• • • 

9.235 

2.8575 

1.27 

8 . . 


• • • 

• 13.335 

2.8575 

1.27 

9 • • 

• • 

• • • 

1.905 

4.7625 

1.905 

10 . • 


• « • 

5.715 

4.7625 

1.905 

11 . • 


• • • 

9.235 

4.7625 

1.905 

12 . • 


• • • 

• 13.335 

4.7625 

1.905 

13 • • 

• • 

• • • 

5-715 

6 .6675 

2.54 

14 • • 

• • 

• • • 

9.235 

6 . 66 75 

2.54 

15 • • 

• • 

• • • 

• 13.335 

6 . 6675 

2.54 

1 6 • • 

• • 

• • • 

. 7.62 

3.81 

0.0 


■ 1 t 


DISCUSSION OF RESULTS 


Three types of fillers \\>ere evaluated in this study 
at five different levels of heat flux. A set of runs 
was made nithout a filler to set a performance base line 
for each power level. The three fillers that were tested 
are presented in table 3> The thickness given in the table 
is the wall thickness of the filler, while the depth is 
the length of the filler in the z-direction. In the theo- 
retical analysis, average physical properties were used 
for n-octadecane, while the literature values, given in 
table 1, were used for the aluminum filler. The average 
physical properties that were used are as follows: 

Density 

= O .8969 gm/cc 

= 0.8545 gm/cc 


p ( solid) 
p ( liquid) 

Heat Capacity 


Cp (solid) 

C (liquid) 

Thermal Conductivity 
K (solid) 

IC (liquid) 


= 2.164 watt-sec/gm/°K 

= 2.406 watt-scc/gm/°K 

= 0.001521 watt/ cm/°I( 

= 0.00735 watt/ cm/°K 


When the literature value of the latent heat of 


fusion was used in the mathematical model, the phase- 
change times were not predicted. This can be seen by 
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Table 3 


GEOMETRY DATA FOR THE ALUMINUM HEXAGONAL FILLERS 


The thickness is the wall thickness of the filler, the cell 
size is the distance across one cell of the filler, while 
the depth is the distance in the z-direction. 


Filler No. 

Thickness 

cm 

Cell Size 
cm 

Depth 

cm 

1 . 

0.00889 

1.905 

2.54 

2. 

0.011938 

0.635 

2.54 

3 . 

0.05969 

0.635 

1.7 
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coinwarin g figures 9 a-d and 25 a-d. Therefore, a reduced 
value of the latent heat of fusion was used. When this 
value was used it predicted- the phase-change times more 
accurately. The latent heat of fusion for n-octadecane 
used in this study is given by Bain (21) and is as follows. 
II„ = 182.83 watt-sec/gm 

If the addition of fillers increases the heat- 
transfer rate into the phase-change material, the hot 
plate temperature should remain below a given control 
temperature for a longer period of time as the amount of 
filler is increased. In figure 10 the ratio of filler 
weight to n-octadecane weight is plotted against a pseudo 
control temperature. This pseudo control temperature is 
the hot plate control temperature minus the initial temp- 
erature, This figure is not intended to be used for 
design, but rather to show that the experimental data 
is consistant . However, if the control temperature minus 
the initial temperature should be 22 K, this figure could 
be used for approximate design purposes. This graph 
should not be extrapolated beyond the experimental data. 

Figures 11 through 34 represent the theoretically 
predicted temperature profiles compared with the 
corresponding experimental data. The thermocouple 
locations are given in table 2. Response data from only 
one thermocouple will be plotted for each height since there 
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Figure 9 

Theoretical and experimental temperature profiles for 
filler no • 2 at 40 watts with the literature value for 
the liquefaction enthalpy 


Figure 

9a 0.635 cm from the heating plate 
9h 1.27 cm from the heating plate 

9c 1.905 cm from the heating plate 

9d 2.54 cm from the heating plate 


Leger. d 

A F2-40-1 

O F 2-40-2 


Theoretical 



ure K Temperature K Tempera 
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is little difference between the temperature profiles 
recorded by each thermocouple for a given height. 

In figures 11 through 16 the theoretical and ex- 
perimental data from the pure paraffin or no filler runs 
are plotted. Due to the low heat-transfer rate of n- 
octadecane and the hot plate temperature limitation 
discussed earlier only a small portion of the n-octa- 
decane melted during these runs. In these figures .the 
paraffin at 0.635 cm has melted and the temperature of the 
liquid is rapidly rising for all power levels. At 1.27 
cm the paraffin has just melted in the 20 and 30 -watt 
runs while it has just reached the melt point in the 40 
watt run. It is still below the melt temperature in the 
50 and 100 watt runs. In the 20 and 30 watt runs the 
paraffin is in the process of melting at 2.905 and 2.54 
cm. In the 40, 50 and 100 watt runs the paraffin at 1.905 
and 2.54 cm is still heating up to the melt temperature. 
The hot plate temperature profiles for the pure paraffin 
runs are presented in figure 16 . 

The filler runs are presented in figures 17 through 
34. The fillers will be discussed in the order in which 
they are presented in table 3> and referred to by the 
number indicated in the table. In figures 17 to 22 the 
theoretical and experimental data for filler number 1 are 
plotted. In these runs more of the n-octadecane has 
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Figure 11 

Theoretical and experimental temperature profiles at 
20 watts for pure n-octadecane 

Figure 

11a 0.635 cm from the heating plate 

lib 1.27 cm form the heating plate 
11c 1.905 cm from the heating plate 

1 1 cl 2.54 cm from the heating plate 


Legend 

A F- 20- 1 
O p-20-2 

Theoretical 


vi *' - r • *i 1 !’‘r {.»• ?• f *• 


» * M > ? j r 


• • 1 


!’IV I 
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Figure 12 

Experimental and theoretical temperature profiles, 
at 3 0 watts for pure n-octa decane 


Figure 

\ 

12a 0.635 cm from the heating plate 

12b 1.27 cm from the heating plate 

12c 1.905 cm from the heating plate 

12d 2.54 cm from the heating plate 

Legend 

A P-30-1 
O P-30-2 


Theoretical 
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Figure 13 

Experimental and theoretical temperature profiles 
at 40 watts for pure n-octadecane 

Figure 

13 a 0.635 cm from the heating plate 

13 h 1.27 cm from the heating plate 

13 c 1.905 cm from the heating plate 

13 d 2.54 cm from the heating plate 

Legend 
A P-40-1 
O P-40-2 
Theoretical 
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Figure 14 

Experimental and theoretical temperature profiles 
at 50 watts for pure n-octaclecane 

Figure 

14 a 0.635 cm from the heating plate 

14 b 1.27 cm from the heating plate 

14 c 1.905 cm from the heating plate 

14 d 2.54 cm from the heating plate 

legend 

A p-50-1 
O r-50-2 


Theoretical 



(c) 
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Figure 15 

Experimental and theoretical temperature profiles 
at 100 watts for pure n-octa decane 

Figure 

15a 0.635 cm from the heating plate 

15 b 1.27 cm from the heating plate 

15 c 1.905 cm from the heating plate 

1 5 <! 2.54 cm from the heating plate 

Legend 

A P-100-1 

O P-100-2 
Theoretical 


( 
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Figure 17 

Experimental and theoretical temperature profiles 
at 20 v/atts for filler no. 1 

Figure 

17 a 0.635 cm from the heating plate 

17 b 1.27 cm from the heating plate 

17 c 1.905 cm from the heating plate 

17 d 2.54 cm from the heating plate 

Legend 

A F 1-20—1 

O Fl-20-2 

Theoretical 



(cl) 



Time sec 





















Figure 18 


Experimental and theoretical temperature profiles 
at 30 watts for filler no. 1 


Figure 

l8a 0.635 cm from the heating plate 
l8b 1.27 cm from the heating plate 
l8c 1.905 cm from the heating plate 
l8d 2.54 cm from the heating plate 

Legend 

A Fl- 30-1 
O F 1 - 30-2 


■ Theoretical 




















Figure 19 
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Experimental and theoretical temperature profiles 
at 40 watts for filler no. d 

Figure 

19 a 0.635 cm from the heating plate 

19b 1.27 cm from the heating plate 

19 c 1.905 cm from the heating plate 

1 9 cl 2.54 cm from the heating plate 

Legend 

A F 1-40-1 

O Fl-40-2 
Theoretical 


P 



Temperature 
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Figure 20 

Experimental and theoretical temperature profiles 
at 50 watts for filler no. 1 

Figure 

20a 0.635 cm from the heating plate 

20b 1.27 cm from the heating plate 

20c 1.905 cm from the heating plate 

20d 2.54 cm from the heating plate 

Legend 
A Fl-50-1 
O Fl-50-2 
Theoretical 
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Figure 21 

Experimental and theoretical temperature profiles 
at 100 watts for filler no. 1 

F igure 

21a 0.6 35 cm from the heating plate 

21b 1.27 cm from the heating plate 

21c 1.905 cm from the heating plate 

2 1 cl 2.54 cm from the heating plate 

Legend 

A FI- 100-1 
O FI- 100-2 

Theoretical 
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melted compared to the pure paraffin runs. The theore- 
tical curves in all cases predict the solid-phase temp- 
erature profiles exactly, while there is some variation 
in the liquid phase. The liquid phase deviations are 
always on the high side. If the theoretical analysis had 
considered the heat losses, the theoretical temperature 
profiles could possibly have been brought down in line, 
with the experimental data. The phase change times are 
all predicted. The theoretical profiles that curve up 
smoothly through the melt point are nodes closer to the 
aluminum filler, figure 17 , while those that jump sharply 
after the phase-change are nodes farther away from the 
filler, figure 17b. The hot plate temperature profiles 
for filler number 1 are shown in figure 22. 

Figures 23 to 28 show the experimental and theo- 
retical results for filler number 2. The same statements 
that were made for filler number 1 can be marie for filler 
number 2. 

The results from filler number 3 are plotted in 
figures 29 to 32. Note that the theoretical temperature 
profiles deviate sharply from the experimental profiles 
at the 1.27 cm level. In this case the filler is only 
1 . 27 cm deep, and when all of the n-octadecar.e has melted 
in the theoretical analysis, the whole cell will heat up 
rapidly as shown in figures 29 b, 29 d, 30 b, 30 d and 31 b. 
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Figure 1 13 

Experimental and theoretical temperature profiles 
at 20 watts for filler no. 2 

Figure 

23 a 0.635 cm from the heating plate 

23 b 1.27 cm from the heating plate 

23 c 1.905 cm from the heating plate 

23 cl 2.54 cm from the heating plate 

Legend 

A F 2- 20-1 

O F 2- 20-2 
Theoretical 
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Figure 24 

Experimental and theoretical temperature profiles 
at 30 watts for filler no. 2 


Figure 

24a 0.635 cm from the heating plate 

24h 1.27 cm from the heating plate 

24c 1.905 cm from the heating' plate 

24d 2.54 cm from the heating plate 


Legend 

^ F 2- 30-1 
O F2-30-2 


Theoretical 
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Figure 25 

Experimental and theoretical temperature profiles 
at 40 watts for filler no . 2 

Figure 

25 a 0.635 cm from the heating plate 

25 b 1.27 cm from the heating plate 

25c 1.905 cm from the heating plate 

25d 2.54 cm from the heating plate 

Legend 

A F 2-30- 1 

O F2-30-2 
Theoretical 
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Figure 26 

Experimental and theoretical temperature profiles 
at 50 watts for filler no. 2 

Figure 

26 a 0.635 cm from the heating plate 

26 b I .27 cm from the heating plate 

26 c 1.905 cm from the heating plate 

26 d 2.54 cm from the heating plate 

Legend 

A F 2- 5 0-1 

O F 2- 50-2 

Theoretical 


e- 
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Figure 27 

Experimental and theoretical temperature profiles 
at 100 watts for filler no. 2 


Figure 



27 a 

0.635 

cm from the heating plate 

27 b 

1.27 ■ 

cm from the heating plate 

27 c 

1.905 

cm from the heating playe 

27 d 

2.54 < 

cm from the heating plate 

Legend 



A F2- 

■100-1 



O F 2-100- 2 


Theoretical 
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A 20 watts 
□ 30 watts 

O 40 watts 
O 50 watts 
© 100 watts 



Figure 28 

Hot Plate Temperature Profiles For Filler Ko . 2 
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This is clue to the fact that the bottom is considered 
insulated and there is no more solid to change phase 
and thus absorb heat. In the experimental cell there is 
definately heat loss through the bottom, which could 
account for seme of this deviation. Also the 2.54 cm 
deep cell was modified to test the 1.27 cm deep filler 
by inserting a 1.27 cm plexiglass plate. It was nece- 
ssary to drill holes in this plate to accomodate the 
existing thermocouples. These holes filled with par- 
affin. With this extra n-octadecane around the cold 
plate or 1.27 cm thermocouples considerable more heat 
is absorbed, thus keeping the thermocouple temperature 
down. The hot plate temperature profiles for the third 
filler are shown in figure 32. It should be noted that 
the inflection in the hot plate temperature profiles are 
probably clue to the presence of air bubbles. 

Figures 33 through 37 demonstrate some of the theo- 
retical studies that can be made with the mathematical 
model developed in this study. In figure 33 the temp- 
erature profiles from various filler wall thicknessess 
are plotted. It can be seen from these plots that as the 
filler thickness increases sc does the heat-transfer rate. 
This is shown by the way in which the slope of the curve 
above the melt point increases as the wall thickness of 
the filler increases. 
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Figure 29 

Experimental and theoretical temperature profiles 
at 20 and 30 watts for filler no. 3 


Figure 

29 a 0.635 cm from the heating plate - 20 watts 

29 b 1.27 cm from the heating plate - 20 watts 

Legend 

A F 3 - 20-1 

O F3-20-2 

Theoretical 

Figure 

29 c 0.635 cm from the heating plate - 30 watts 

29 d 1.27 cm from the heating plate - 30 watts 

Legend 

A F3-30-1 

O F3-30-2 


Theoretical 
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Figure 30 

Experimental and theoretical temperature profiles 
at40 and 50 watts for filler no. 3 

Figure 

30 a 0.635 cm from the heating plate - 40 watts 

30 b 1.27 cm from the heating plate - 40 watts 

Legend 

A F 3 - 40-1 

O F3-40-2 

Theoretical 

Figure 

30 c ' 0.635 cm fro the heating plate - 50 watts 
30 d 1.27 cm from the heating plate - 50 watts 
Legend 
A F3-50-1 
O F3-50-1 


Theoretical 











85 


Figure 3 i 

Experimental and theoretical temperature profiles 
at 100 watts for filler no. 3 

Figure 

30 a 0.635 cm from the heating plate 

30 b 1.27 cm from the heating plate 

Legend 

A F3-100-1 
O F3- 100-2 
Theoretical 
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/S 20 watts 
□ 30 watts 
O 40 watts 
O 50 watts 


(J) 100 watts 



0 1200 2400 3600 


Time sec. 

Figure 32 

Hot Plate Temperature Profiles for Filler No. 3 
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Figures 34 through 37 are included to show the 
temperature profiles in the x-y plane. In figure 34> 
there is not much deviation. in the x-y plane. This is 
intuitively correct since the graphs were made from data 
taken at a time early in the run. As time increased so 
does the temperature variation in the x-y plane. (Fig- 
ures 35, 36 and 37) Kote that the n-octadecane closer 
to the filler rises to a higher temperature than that 
farther away at any given time. This type of heating 
profile is exactly what is expected with a high thermally 
conductive metal matrix in a phase-change environment. 

The experimental data in this section could be used 
for design providing that the design requirements fall 
within the experimental data and situation as presented 
in this study. The computer program written for this 
study can be used to predict the capabilities of other 
thermal-control devices by varying the physical prop- 
erties of either the filler or phase-change material. 

If a different filler geometry is to be studied the 
subroutine which calculates the areas of the filler and 
phase-change material must be changed to accommodate the 
different geometry. If heat loss from the system is to 
be considered, the computer program can be modified 
as described in the theory section to provide for this. 
The computer program written for this study can also be 
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Theoretical Temperature Profiles for 
Different filler Thicknessess 
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Figure 34 


Theoretical temperature profiles in the x-y plane at t= 
600 sec 

c 
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Figure 36 

'V . 

K'. m . 

Theoretical temperature profiles in the x— y plane at t= 

1800 sec 
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Figure 37 


Theoretical temperature profiles in the x-y plane at t= 
2400 sec 
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used to predict the performance of thermal-control 
phase-change devices under different heat loads by 
modifing it to use a heat-flux boundary condition a 
described in the theory section. 



RECOMMENDATIONS AND CONCLUSIONS 


Conclusions 

Based upon the results of this study, the following 
conclusions are presented. 

1. As the weight of the filler material increases 
the heat-transf er rate of the thermal-control phase-change 
device increases . 

2. The computer program written for this study 
predicts the experimental solid-phase-temperature profiles 
and the phase-change times correctly. While the maximum 
deviation between the theoretical and experimental solid- 
phase-temperature profiles is 2.8°K (5°F), there is 
essentially no deviation in most of the runs. There is 

as much as 14°K (25°F) deviation between the theoretical 
and experimental results in liquid phase. This deviation 
could possibly be corrected by changing the boundary 
condition along the bottom plate, in the computer program, 
to one in which heat loss is allowed. 

3 . The mathematical model presented in this study 
is general in terms of variable filler geometry, physical 
properties of the filler and phase-change material and the 
types of boundary conditions that can be placed on the 
theoretical model. The computer program written for this 
study is general in terms of variable physical properties 
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of the filler and phase-change material and sizes and 
shapes of the hexagonal filler. While the computer 
program is written for insulated boundaries along the 
bottom and four sides, it can be easily modified to 

\ 

allow for heat loss or gain along these boundaries. 

The computer program uses a temperature profile on the 
heating plate but this boundary, as discussed earlier, 
can also be changed to utilize a heat-flux boundary 
condition . 

4. A three-dimensional analysis is needed’ to see 
the detailed temperature profiles in the x-y plane. As 
shown in the discussion of results section these theo- 
retical temperature gradients can be as much as 11.1°K 

( 20°F ) . 

5. This study has shown that to correctly model 
the experimental system the heat losses must be known 
or predicted. If a heat-flux boundary condition is to 
be used in the theoretical model, the heat-flux into the 
test chamber must be accurately known. 

Recommendations 

The following recommendations are presented based 
on the results of this study. 

1 . The test chamber should be redesigned in such 
a way as to eliminate air bubbles. 

2. Since the n-octadecane used in this study tended 



to trap air in the solidification process, other phase- 
change materials should be considered. These other 
materials could include lithium nitrate trihydrate and 
acetamide . 

3. Since the filler material adds weight to the 
phase-change thermal-control unit, and subtracts from 
its heat-absorbing capacity, other filler materials and 
geometries should be studied with the goal of optimizing 
the ratio of filler material to phase- change material. 
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NOMENCLATURE 


a,b, c 


A 


C 


P 


G 


h 


H 


f 


I 


J 


K, 


K 

q 

T 


Coefficient of the unknown temperatures in the 
tridiagonal matrix 

The cross sectional area perpendicular to the 
2 

heat flux; cm 

The heat capacity; watt-sec/gm/°K 

Amount of energy generated per unit volume; 
watt-sec/ cc 

2 

The heat-transfer coefficient; watt/sec/cru 

Enthalpy of .liquefaction; watt-sec/gm 

The last node in the x-direction in the nodal 
network 

The last node in the y-clirection in the nodal 
network 

The last node in the z-direction in the nodal 
network 

Thermal conductivity; watt/cm/°K 

2 

The heat flux; watt-sec/ cm“ 

Temperature; °K 
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T 

e 

Th 


T 

mo 

T 

m 

t 

At 


V' 


V 


P 

<P C p V) eff 


(^ A > ef f 


The excess degrees; °K 
Wall thickness of the filler; cm 
The initial melting temperature:; °K 
The melting temperature; °K 
Time; sec 

Incremental time; sec 

The volume of the material that generates 
energy; cc 

Volume; cc 

Dens ity ; gm/ cc 

The effective (pC^V) for a nonhemogeneous node; 
watt- sec/°K 

The effective (KA) for a nonhomogeneous node 
watt-cm/°K 


Subscripts 

x X-direction which is along an axis parallel to 

the 15.24 cm side of the test chamber 

y Y-directiors. which is along an axis parallel to 

the 7.62 cm side of the test chamber 
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z Z-direction which is along an axis parallel to 

the 2.54 cm side of the test chamber 

i Indicates the x-direction in the finite 

difference formulation 

j Indicates the y-direction in the finite 

differeiice formulation 

k Indicates the z-direction in the finite 


difference formulation 



APPENDIX I 


Computer program to solve the nonhomogeneous phase-change 
problem. 

This computer program was written in FORTRAN IV to 
solve the nonhomogeneous phase-change problem presented 
in this study. This prograin uses 14 cards of input data. 
The coefficients of the three straight-line fits for the 
hot plate temperature profiles are read on the first 2 
cards. The physical properties of the filler and phase- 
change material are read on the next 9 cards. The filler 
geometry is specified by the last 2 cards. The following 
table specifies the exact variable to be read on each card. 


VARIABLE 
Card 1 

A1 Coefficient of t in T = At + B; °F/min 

B1 B in T = At + B; °F 

Til The last time for which A1 and B1 hold; min 

Card 2 

A 2 Same as A1 

B2 Same as B1 

TI2 The last time for which A2 and B2 hold; min 
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VARIABLE 
Card 3 
A3 
B 3 

Card 4 
DIF 
TM 
TA 

Card 5 
CPM 

HF 

TS 

Card 6 
ATI 
XIC 
Card 7 
TH 
RF 
CPF 

Card 8 
FK 

Card 9 
DX 
DY 
DZ 


Same as A1 
Same as B1 

The melting temperature range; °F 
The initial melting point; °F 

o 

The initial temperature of the test cell; F 

The heat capacity of the melting paraffin; 
BTU/lb/F 

Enthalpy of liquefaction; BTU/lL 
The ending time of the run; min 

The time for the first print cut; min 
The time increment between print outs; min 

The wall thickness of the filler; in 

Density of the filler; lb/cf 

Heat capacity of the filler; BTU/lb/°F 

Thermal conductivity of the filler; BTU/ft/hr/°F 

Length of the node in the x- direction; in 

Length of the node in the y-direction; in 

Length of the node in the z- direction; in 
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VARIABLE 
Card 10 
DT 
PKS 

PKL 

Card 11 
RPS 
RPL 

CPS 

Card 12 
CPL 

Card 13 
NX 

11 

12 

Card 14 
NY 
N2 
N4 


Time increment; sec 

Thermal conductivity of the solid phase-change 
material; BTU/ft/hr/°F 

Thermal conductivity of the liquid phase-change 
material; BTU/ft/hr/°F 

Density of the solid phase-change material; lb/cf 

Density of the liquid phase-change material; lb/ cf 

Heat capacity of the solid phase-change material; 
BTU/lb/ F 

Heat caoacity of the liquid phase-change material; 
BTU/lb/°F 

The number of nodes in the x-direction plus 1 

The first node of the angular section of the filler 

The last node of the angular section of the filler 

The last node in the y- direction plus 1 
Set equal to 3 

The last node in the z-direction plus 2 
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DIMENSION VP<7,7> ,VF<7,7> 

D I HEMS I ON A(32> ,n(32) ,CC(32) |D(32 ) iT(32) 

DIMENSION A P X ! ( 7 , 7 ) , A f ' X 1 ( 7 » 7 ) 

DIMENSION APY! (7, 7) ,AFYI (7,7) 

DIMENSION T1(7,7i3?>iT2(7,7,32> i T3(7i7,32> 

DIMENSION TE<7, 7, 32), KJ(7, 7, 32) 

DIMENSION APH ( 7 > 7 ) » AF? ( 7 , 7 ) 

COMMON /SI/ N4, I , JiK 
COMMON /$?./ VP,VF 
COMMON /S3/ APXI,AFXl,NX 
COMMON / S 4 / A P Y I , A P Y I i N Y 
COMMON /S3/ A P / , A F 2 | N 2 
COMMON /S6/ A , B , CC i 0 
COMMON /S7/ TM l HF|CPM l DIF f VT 
COMMON /SB/ Al, A2|A3 iB1iB 2»B3,TU,TI2»TI 
COMMON /S9/ PKS,PKU 
COMMON /si'// RPL,« p S,CPS,CPL 
,|,,,,N2 NODE THAT STARTS THE TEST M A TER I AU 
• » • i » t N 4 * the no, OF NODES in THE TEST MATERIAL 
, . i i . i VF= VOL OF FILLER * ROE ( F 5 *CP ( F ) /DT 
M.,,. VPsVOL OP PARA, » ROC(P)»CP(P)/DT 
, , . , , ,NXs NO, of nodes in x dir, 

, . , , , ,NY= NO, op nodes in y oi r, 

. , , , , ,MEr NO, OP NODES JN ? D I R , 

I , . , Vi TM :• MELT TEMP. OF PARA, 

. , . , . ,MFb' 'heat OF FUSS I ON OF THE PARA, 

. . . cpnshfat capacity or the para, at the melt temp, 

. . . , . .CPPeUEAT CAPACITY of para, 

. , . , . .CPFeUEAT CAPACITY OF FILLER 
, , • , ', » RFs DENSITY OF FILLER 
. . , , p P = DENSlTY OF PARA, 

, , , , , .PKsTHERMAL COND, OF PARA, 

. , , . , .FKcTHERMAL COND, OF FILLER 
READ (2. 6/2) A1 , Dl , T 1 1, A2 , B2 , T 1 2 , A3 , B3 
READ (2. 600) DIF,TM#TA,CPM,HF,TS,ATl,XlC 
R E A Q ( 2 , 6 0 0 ) TH,RF,CPF,FK 
READ (2, 603) DX , DY , D2 , DT , PKS , PKL 
READ ( 2 » 6O0 ) RPS, RPL, CPS, C p L 
READ ( 2 i 60 1 ) MX, II, 12»MY,N2,N4 
609 FORMAT (61) 

60P FORMAT ( 3F ) 

601 FORMAT (31) 

WRITE (1,602) Al,!31»TIli A2.B2, TI2, A3,B3 
WRJTE(1,603) DIF,TM,TA,CPM,HF,TS,ATI,XIC 
WRITE(1,604) TII,RFiCPF,FK 
WR!TE(1,605> DX,nY,Q/,DT ,PKS,PKL 
WRITE (1,607) RPS,. RPL, CMS, C p L 
WRITE (1,606) MX, II, I2,MY,N2,N4 

602 F0RMAT(4X, 'Al=* ,P0,4, t DEGF/MJN' ,7X, *Bls» ,F8.4, • DEGF ' t 
15X, «TU = » ,F0, 4, » MJN*/,4X, * A2=' »F6,4, ' DEGF/MIN* 

27X, »B2=' ,Fl,4, ' DEGF'.SX, 'TI2 S ',FC.4, » MIN'/, 

34X, * A 3 = * ,m,4, » D E 0 F / M I N 1 , 7 X , ' D 3 = ' , F 8 , .4 , ’ DEGF') 

603 FORMAT ( 3X , * D I F= 1 , F8 , 4 , ’ DEGF 11X , ' TM= * , F6 , 4 , » DEGF’,6X,'T 
IF 8*4, ’ OFGF’ ,/3X, ♦ C P M = « , F 3 , 4 , * B T U / L B ' ,9X, ' H F = ' , F 8 , 4 , ' BTU/L 
24X, »TS=» ,FB,4, » MI'JI ,/,3X, 'ATI = ' ,F8,4, ' MIN' , 1 1 X , ' X I C = * , 

3F8.4,' MIN') 

604 F0RMAT(4X, , TH=',F0.4, » I NCHES « , 9X , ' RF= » , F 8 V 4 , ’ LB/CF ' , 4X , ' CPR= ' 
1F8.4, ' BTU/LB/DEGF' ,/,4X» ’ F K = ' ,F8,4» ' BT'J/FT/MR/DEGF ' ) 

605 FORMATUX, «DX=« ,F8,4, » INCHES », 9X» »0Y® ' ,F8, 4 , i INCHES’ 

1#4X, 'DEa' ,FB,4, ' I NCHES » , / , 4 X , ’ 0T = • i FO , 4 , • SEC , » , 4X , ’ PKS = * , F 3 , 4 , 


C, , , 
C i i , 

c , » , 

C i • , 

C , i I 

c , , , 
c , , , 
C i . i 

C , i i 
C , , i 
C , , a 

C I < a 

Ca, a 
Ci 1 1 
C i i a 


< m 
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1«8TU/FT/HR/0EGF« ,3X, ' PKL = * , F3 t 4 , * BTU/FT/HR/DEGF » > 

606 FORMAT MX, 1 N X s * i I 3 » 2 1 X # ’ II 1 ? 1 » I 3 , 1 6 X , • I 2 *= * , I 3 # / 4 X # * N Y = 1 i 13# 

1 2 1 X i *N2?» , 13»16X» «M4»' , 13) 

607 FORMAT (4X, ' RPS= ' , FQ , 4 , » L8/CF * , 3X , • RPL = ’ FB , 4 , ' LB/CF ' 

13X, * CPS?? « ,F8,4, ' BTU/UB»,/»4X, 'CPL=» ,F0,4, 1 BTU/LB') 

599 N3«»N2*1 

NZH"N4*1 
TI*0, 

N2 = N4 

DT=OT/60,/60, 

NXbMX+1 

NY?MY*1 

DX«DX/12,' 

DY=DY/12 , 

02=02/12, 

VT = QX«>DY«QZ 
TH=TH/12, 

CAUL AREA( II. 12.TH»DX,DY,DZ, DT,RF»CPF»FK> 

DO 5 K = 2 » N 2 
DO 5 J = 2 , N Y 
DO 5 lag.NX 
T 1 ( I , ,J , K ) = T A 
T2U,J,K)sTA 
T3(!iJ,K)=TA 
KJ(I,J,’K>«1 

5 CONTINUE 

6 CALL TBN(TBl) 

TJ *TI*2 , *OT*60V 
CALL TBNCTB2) 

TB=<TBl+TB2>/2, 

DO 7 J = 2 , N Y 

DO 7 1=2, NX 
TiU, J,2)#TB 
T2 ( I , J , 2 ) =TB 
T3 U , J , 2 ) = Tg 

7 CONTINUE 
CVVV'.i.t SOLVING FOR T2 

DO 15 K = 3 1 N Z 
KK=K+1 

DO 15 J = 2 , N Y 
JJ-J+1 

DO 10 1=2, NX 

CALL KC(C,KJ{ I, J, K > > 

11=1+1 

CALL KKCll# I ,KJ(I , J,K> ) 

CALL KKC12, 1 1 ,KJ< I , J,K) ) 

CALL K2(C21, J,KJ( I, J,K> > 

CALL K2 ( C22 i J J i K J ( I , J , K ) ) 

CALL K?(C31,K,KJ( Jf J,K) ) 

CALL KS ( C32 , KK i K J < I , J » K ) ) 

621 FORMAT! » CS»,3I2,4E) 

610 FORM AT ( 1 COE ’ 3 1 2 , 5E ) 

AU >*sCll/C/DX 
BU )=l, + (Cli + C12)/C/DX 
CC(I >s-cl2/C/DX 
D ( I > = T1 ( I , J,K ) 

ZlP(C2i»Tl(Ii JaliKI-(C21*C22)*Tl( I, J,K)+C22«T1< I. J+1,K) )/C/DY 
D(I)=D(I)+Z1 

Z2P(C31*T1< I, J,k;i>»(C 31*C32)*T1< l, J,K)*C32*T1( I, J,K*i) >/C/DH 
612 FORMAT ( 3X , 4£ > 


-i ■ ■ « ■ i-i n v, « irrri , 


(' ■,*»* - ,'n it •[ t , . »• i ~ %. r i * *, < i I r »• - i *■,«, 


1 ] 2 

DU>=D(I>*22 
611 FORMAT ( 3E ) 

10 CONTI NUF 

CALL SOLVE < T , 2 , NX > 

DO 12 K I a2 » NX 
12 T2 ( K I , J » K ) =T ( K I ) 

15 CONTINUE 

C| , Vi t i SOLVING FOR T3 
DO 25 I “2 i NX 
11=1+1 

DO 25 K b 3 1 N 2 
KK*K+1 

DO 20 J*=2»NY 

CAUL KC(C.KJ(I,J,K>) 

JJ = J + l 

CALL K2(C21.J»KJ( I, J»K> ) 

CALL K?(C22#JJ»KJ( I. J,K> ) 

A { J > - C21/C/DY 

B( J)b ,.*<C21+C22>/C/DY 

CC( J>=rC22/C/DY 

H«(C21<-C22)*T1(I, J»K).C21*T1( I, J«1 , K > «C22*T1 ( J,J*1,K> 
ZpH/O/DY 

D(J)"T2( IiJ»K>*Z 

20 CONTINUE 

CALL SOLVE (T,2,NY> 

DO 22 K I a2# NY 
T3(I,KI,K)=T(KI) 

22 CONTINUE 

25 CONTINUE 

DO 28 Ks 2,NH 
DO 28 U=2»NY 
DO 28 I ?2 » NX 
T2n»J,K)sTl(JiJ,K) 

28 CONTINUE 
cVVVVi • .SOLVING FOR T 1 

DO 35 J=2jNY 
JJ=J+1 

DO 35 1=2, NX 
11=1+1 

DO 32 K * 3 1 N 2 

CALL KC < C » K J ( I , J , K ) ) 

KK=K+1 

CALL K3(C31,K,KJ( ) 

CALL K3(C32,KK,KJ( l , J.K> ) 

C- 3 3 = C 3 1 
EC = 0 , 

I F < K , GT . 3 > GO TO 29 
ECsT2( I , J#K-1)«C31/C/DZ 

29 A(K>Uc33/C/DZ 
8(K>=1,+(C31+C32)/C/DE 
CC<K)s s C32/C/D?. 

D(K>52,#T3(I,J,K)wT2(I,J,K> 

H?{C3l*C32>*T2(I,J,K>»C3l*T2n,JiKii)^C32*T2(IiJ,K+l> 

Z*?./C/DZ 

D < K ) =D(K ) + Z*EC 

30 CONTINUE 

CALL SOLVE ( T , 3 * NZ ) 

DO 32 KI*3|NZ 
Tl( I# J,KI)=T(KJ> 


H3 


32 

35 


37 

38 


40 
444 
45 

41 

42 

43 


999 


CONTINUE 
CONTINUE 
* DO 37 KsM2iNH 
DO 37 J=2»NY 
DO 37 I 5 2 NX 

CALL PHASE? Tl ? I , J, K) , KU? 1 1 J»K) , TE ? I # Ji K > , VP? 1 1 i T2 ( I , J, K > > DT ) 

CONTINUE 

IF C T I . U T . A T I ) GO TO 6 
WR I TE ( 1 i 42 ) TIiTQ 
DO 40 K=2,NH2,4 
L1 = K<-1 
L2=K+l 

WRITE?-: 41) LI, (Tl? I,2,Ll) » 1*2. NX) 

WR I TE ( 3 , 41 ) Ll , ? Tl? I * 3 , Ll ) i 1 = 2, NY) 

WRITE?!, 4l> K,?Tl(l*2»K),I=2,NX> 

WRITE?! .41) K , ? Tl ( I » 3 , K ) , I = 2 , NX ) 

WRITE (1,41) L2,?T1?I,2,L2),I*2,NX) 

WRITE? 1,41) L2, ?T1< I »3,L2) » I =2, NX 5 
WRITE (1,444) 

CONTINUE 

format?//) 

F0RNAT(4(3X,3I2,F6,1) ) 

FORMAT(3X, I3,4(3X,F6,1> ) 4 

FQRMAT(//3X, ' T I H E 8 ’ #E8,3, ’ MIN,’,3X,'TB= , ,E6 I 1*/) 

FORMATdlX, 12(212, 5X)) 

ATI-ATI+XIC 

IF(TI.GT.TS) GO TO 999 

GO TO 6 

STOP 

END 

SUBROUTINE AREA ('Hi 12, TH,DX, DY, DZ, DT,RF *CPF »FK) 

DIMENSION VP (7, 7) i VF(7»7) 

DIMENSION APXI (7,7) , AFXI (7,7) 

DIMENSION APYI (7,7) » AFYI ?7i75 
DIMENSION Af 5 2(7,7)iAF2?7»7) 

COMMON /S2/ V P , V F 
COMMON /S3/ APXJ , AFXI ,NX 
COMMON / S 4 / APYI,AFYI,NY 


COMMON /S5/ APH,AFZ,NZ 
, , ,DX>OY/(TAN?AA«90,DEG I ) 

, V V. .7. WHERE] 

,,,,,,,, I, AA5 THE ANGLE THE FILLER MAKES WITH THE H 0 R | S £ 0 N T A L 
, ,V. iII-NO, OF NODES IN THE X DIR, 

,'JJaNO, OF NODES IN The y DIR, 

.......... 11= MO, OF NODES TO THE ANGLE PART OF THE FILLER 

I?= NO. OF NODES IN THE ANGLE PART OF THE FILLER 
,,.,,.,...11 5 12 MUST BE INTERCERS 
, ,V, , , V.THsFILLER thickness 

,OX=DELTA X 
.I.OYpDELTA Y 
.V.DHbOELTA Z 

THIS SUBROUTINE WILL PRODUCE THE FOLLOWING ARAYS] 

. .V, ..... ,’APXIsAREA OF PARA, in node 1 1 J IN THE X D I R » 

AFX Is APE A OF FILLER IN NODE I,J IN THE X DIR, 

,apy i parea of par a, in node i,j in the y dir, 

. .".AFYIbAREA OF FILLER IN NODE I , J IN THE Y DIR, 

; , , , .APZsAREA OF PARA IN THE H D I R « IN NODE I , J 

iAFZp AREA OF FILLER IN THE H DIR, IN NODE IiJ 

vp= vol. of para in node i,j 

VF?VOL. OF FJLLER IN NODE I , J 
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IUsNX*i 
JUaNYoi 
I L c 1 
JL*1 

TAflFLOAT< 1 2-11) /FLOAT (JU.JL + l) 

A A ? A T A N ( T A > 

RsFLO AT ( I2j.I 1 )*0X*FLOATUi>*DX' 

00 53 J=JL,JU 

YL?rL0AT( jU„J*i) tt QY 

YU?FU0AT( JU*J)«DY 
XL = S«TAttYl» 

XUb9<?T A »YU 
DO 53 1 ~ I L # I U 
IFCI.Lt.Xl> GO TO 10 
IF< I ,LE, 12) GO TO 20 
! F < I , Lt . I U > GO TO 33 

10 X F < J , GT , JL > GO TO 14 

11 APXI(I,J)?DE«(DY«TH/2,> 

APY I ( I » J ) s0 , 

AFE( I , J)«0X*TH/2, 

AFX I < I, J)bD2*TH/2, 

AFYI(I,J>?0, 

IFd.EO.IL> GO TO 12 
!F< J»EG, JU) GO TO 34 
GO TO 52 

12 APXKI,J)s0, 

AFX I C I , J ) ?0 , 

GO TO 52 

14 APXI ( I , J)=DY«Q2 

APYI ( I , J)*0X*02 
A F X I ( I , J ) « O , 

AFY I ( I , J > s0 , 

AF2 ( I f J ) s0 , 

IF ( XBL . tQ , XL > GO TO is 

if<xbl*dx,e:q,xu> co to 19 
16 IFU.EO.IL) APXI<I,J)=0, 

IF(J,EO,JL) APYX < I , J>»0, 

GO TO 52 

18 APYI(I,J)s 02#(DX'»TH/4 i > 

AFY I ( I # J ) 3DX«TH/4 1 

GO TO 16 

19 APXI ( I , J)sD2«(0Y;TH/4, ) 

AFXJd, J)b02*TH/4, 

GO TO 16 

20 XBLsFLOAT ( I )<*DX 
IF(X»L,LE,XL) GO to 14 
IF(XBUDX,GE,XU> go TO 14 
XlBXBUnXL 

I F < X 1 , G T , D X > GO TO 28 
IFtXU.LT.XBL) XlsXUwXL 
H18X1 /SIN( Aa) 

yi«xi/ta 

APXI ( I , J)=D2#DY 
APYI ( I, J)sD2*(0X;TH) 

AFZ ( I , J) sHl»TM 
AFX I ( I, J)s0, 

AFY I ( I, J)an2*TH 

21 I F < XBL^DX , EQ , XL ) GO TO 31 

22 IFd-lVE'Ml.ANDVJiEO.JL) GO TO 26 

23 IF < J, EO , JL ) GO TO 25 
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GO TO b 2 

25 ARY I ( I , J)eci, 

AFYI ( 1 1 J)3fi, 

GO TO 52 

26 APXI ( I , J>bD2«(DY;TH/?, ) 

A T X 1 ( 1 , J)cOZ*TH/;>, 

GO TO 23 

28 Y2pDYhY1 

H1pY2/COS ( A A > 

A PX I ( I,J>=DH*(OYfTH) 

APY I ( I , J) = OXOZ 
AFZ( I, j) = TH«ni 
AFX I ( 1 , J>bOZ*TH 
A F Y I ( I , J ) = O , 

GO TO 21 

31 APXKl, J)sOH*(nY-TH/4,> 

AFX I ( I, J)sOF*TH/4, 

APYI ( I , J ) = 0 2 « ( n x 1 T H / 4 , > 
AFYI (I, J>sAFXf(I,j) 

GO TO 22 

33 IFCJ.ECi, JU> GO TO 11 
GO TO 14 

34 aPYJ ( J ,J) =OZ*DX 
GO TO 52 

52 APZ(I, j)sDX*PY*AFH< I, J) 

VP ( 1 1 J ) =/.P2 ( I , J ) 02 
VF(Iiv))bAFH(I, J)«02 
APXI (IU+1, J)=2, 

APYI ( I , JU*l)s2l, 

AFX I ( I !.;♦!, J):?, 

APYI ( I , JU*l)s0, 

XBU=0, 

XIO, 

VF ( I, J)sVF(Ii J)OF»CPF/DT 
VP ( I , J)sVP( I, J)/ dT 
AFX I ( I, JJsAFXl < I , J)#FK 
AFYI ( I # J ) s A r Y I { I , J ) <t F K 
AFZ(I,J>sAF2(I,J)ttFK 

50 CONT ! IJUF 

100 FORMAT ( 1 X < 2 1 2 1 4 £ ) 

101 FORMAT C IX, 4£) 
APXI(IU*l»JU*l)s?, 

APYI ( IU*1, JIJ+i ) s2, 

AFX I ( IU*1, JU*1) = M, 

AFYI { I U + l , viU*l ) = 3 , 

DO 200 I = I U , It.,-1 
DO 200 J=JU,JL,-1 
APXI(I*1,J*1)BAPXI(I # J) 
APYI < I*l» J*l>sAPYl ( I , J) 
AFXI ( i*i, J*1 )saFvI ( I , J) 
AFYI ( 1*1, JH> = AFVI < I « J ) 
APH{ I + i, j+3 ) = Af’2 ( I i J > 

AFZ( 1*1, J*l)sArz ( I * J) 

VP( 1*1, J*1>BVP( I , J) 

VF(I* 1 , J+l)sVF(I, J) 

200 CQNTINl'L' 

DO 210 1 = 1, I U 
APX I ( I , JL ) =0 , 

AFX I ( I , JL ) =■” , 

APYKI, JUJPtr, 
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AFYI ( I, JL>*0, 

AP2( I, JL)p0, 

AFZ( I, JL)b0, 

V P ( I i J L ) a 3 , 

VFU, Jl)s0» 

213 CONTINUE 

00 215 Jal.jU 
APXI(IL,J>?0 l 
AFXJ < It, J)a0, 

APYKIU. J>*0, 

AFYI ( Iti J)=0, 

APZdLi J)a0, 

AFH( IU, J>o0, 

VP< JL, J)b0, 

VF( II, J)e0, 

215 CONTINUE 

501 F0PHAT(2X, 212, 2F) 

RETURN 

END 

SUBROUTINE KC ( C i K2 ) 

DIMENSION VP (7, 7) i VF < 7 # 7 ) 
COMMON /SI/ N4 #I,JiK 
COMMON /S2/ VP , VF 
CAUL HC<K2,RP,CPP) 
vi«RP»r.pp*vp( i , j) 

V2? VF ( I # J ) 

C b V 1 + V 2 

I F ( K , E 0 , N 4 ) CsC/2, 

100 FORMAT ( 1 SC*,3E> 

101 FORMAT ( * SC 1 , 3 J 2 , 4E ) 

RETURN 

END 

SUBROUTINE K1 < Cl , 1 1 , L > 
DIMENSION APXI (7,7 ) » a F X I ( 7 1 7 ) 
COMMON /SI/ N4, I , J»K 
COMMON /S3/ apxi ( aFxi,nx 
IF ( I r*E0,N v ' + l> GO TO 63 
I F ( I I . f Q . 2 ) GO TO 60 
CALL T B ( L » P K ) 

A 1 ? A P X I ( I I , J ) #PK 
A2 = AFX I ( 1 1 » J) 

I F ( K , E Q i N 4 ) GO TO 40 
C V V ,, .PARAFFIN AND FILLER NODES 
GO TO 65 
40 AlsAl/2, 

A2=A2/2, 

GO TO 65 
60 A 1 a 0 , 

A2®3 i 

65 CiaAl*A2 

100 FORMATC SCI’, 312. 2E) 

101 FORMAT (4E) 

RETURN 

end 

SUBROUTINE K2 { C2 . JU » L ) 
DIMENSION APY I ( 7 » 7 ) » AF Y I { 7 , 7 ) 
COMMON /SI/ N4, I * J , K 
COMMON /S4/ A P Y I , A F Y I , N Y 
I F ( JJ.EQ.NY+l) GO TO 60 
I F ( JJ , EQ , 2 ) GO TO 60 


c-- 
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CALL TML.PK) 

A 1 3 A P Y I ( I , JJ)*PK 
A2=AFY! ( I » JJ> 

I F ( K , E 0 , 4 ) GO TO 40 
CV'.iV.'. PARAFFIN AND FILLER NODES 
GO TO 65 
40 A 1 a A 1 / ? , 

A2a A2/2 i 
GO TO 65 

60 Ai«0, 

A2a0, 

65 C2» Ai* A2 

100 FORMAT ( ’ SC2',3I2,2E) 

101 FORMAT ( 4E ) 

RETURN 

END 

SUBROUTINE K3 ( C3 1 KK * L ) 

DIMENSION AP2(7»7> t AF2(7i7) 

COMMON /Si/ M4i I . JiK 
COMMON /S5/ APHiAFHiNH 
EC»3. 

I F ( K K , e; Q , N 2 * 1 ) GO TO 60 
CALL TH(LiPK) 

A 1 b A P 2 ( I , J)*PK 
A2»AF?U,J> 

CVVV'.V .PARAFFIN AND FILLER NODES 
GO TO 65 
60 A 1 a 3 * 

A2a0 , 

65 C3=A1*A2 

100 FORMATC SC3’,3I2»3E> 

101 FORMAT ME > 

RETURN 

END 

SUBROUTINE SOLVE (T,J1,J4> 

DIMENSION A(32) »B(32) »CC(32) ,0(32) »T(32) 
DIMENSION W ( 3 2 ) , G ( 3 2 > 

COMMON /S6/ A , B , CC i 0 
J2 S Ji+1 
U3=J4«1 
W(J1)bB(J1> 

G( Jl)sO( Jl)/W( Jl> 

DO 23 I=U2,U4 
XsCC< It.l>/W(I.i) 

W ( I ) =B ( I > a A ( I ) «X 
G( I >s(D( 1 >*A< I )#G< Nl) )/W( I 5 
20 CONTINUE * 

T(J4)=G(U4) • 

DO 33 IaJ3,Jl,si 
XaCCU >A'< J ) 

30 T( I >=G< I >«X*T< I +i > 

101 FORMAT (/> 

100 FORMAT (5E> 

return 

END 

SUBROUTINE PHASE(T,KZ,TX,V,TO,DT) 

COMMON /'67/ TM , HF * CPM | D I F » VT 
VPbV^OT 

J F ( K 2 , G T » 1 ) GO TO 20 
IF < T , LT , TM > GO TO 20 
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SUM = TtsTO*TX 

. ir(SUH*VT»CPM/VP t CE l HF) GO TO 12 
T X ? SUM 

T*TM*DJF«SUM#CPM»VT/VP/HF 
GO TO 20 

12 T?DIF*TM*(SUH*CPM*VT/VP*»HF>/CPM 

K5f»2 

20 RETURN 

END 

SUBROUTINE TH(LiRK) 

COMMON / 59/ PKS.PKU 
IF(L,GT,1) go to 10 
PKpPKS 
GO TO 20 
10 PKpPKL 

20 • CONTINUE 

RETURN 
END 

SUBROUTINE HC(t»RP»CPP> 

COMMON / S 1 O / RPt,RPS,CPS,CPL 
I F ( L i GT , 1 > GO TO 12 
CPPsCPS 
RPsRPS 
GO TO 14 
12 RPsRPL 

CPP-CPL 
14 RETURN 

END 

SUBROUTINE TBN(TB) 

COMMON /SS/ A1,A2* A3»81,B2»B3iTU,TI 2»TI 
IF < T I .lE.TU) GO TO 10 
IF ( T I V L E V T 1 2 ) GO TO 20 
GO TO ?0 

10 T8«A1*TI*B1 

GO TO 40 

20 TB = A2<»T I + B2 

GO TO 40 

30 TB?A3#TI*B3 

40 CONTINUE 

RETURN 
END 



